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We discuss realization, properties and performance of the adaptive finite element approach to the design of 
nano-photonic components. Central issues are the construction of vectorial finite elements and the embed- 
ding of bounded components into the unbounded and possibly heterogeneous exterior. We apply the finite 
element method to the optimization of the design of a hollow core photonic crystal fiber. Thereby we look 
at the convergence of the method and discuss automatic and adaptive grid refinement and the performance 
of higher order elements. 
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1 Introduction The growing complexity and miniaturization of nano-optical components makes ex- 
tensive simulations indispensable. Since in modern devices and applications the wavelength of light is of 
the same order as the dimension of the simulated structures Maxwell's equations have to be solved rigor- 
ously in order to get accurate results for the electromagnetic field. Examples of such structures are meta 
materials, photonic crystal devices, photolithographic masks and nano-resonators lO |6l |7] [lOl (TS) . A lot 
of different simulation techniques have been applied to and developed for nano-optical simulation, e.g. the 
finite element method (FEM), finite difference time domain simulations (FDTD), wavelet methods, finite 
integration technique (FIT), rigorously coupled wave analysis (RCWA), plane wave expansion methods 
(PWE). 

Here we present the finite element method for the solution of time harmonic Maxwell's equations, i.e. 
we compute steady state solutions for electromagnetic fields with a fixed frequency uj. 

This article is structured as follows. First we look at different problem classes which correspond to 
typical nano-optical simulation tasks and give the corresponding mathematical formulations of Maxwell's 
equations. These are propagation mode problems (Sec. [3]), resonance problems (Sec. HI and scattering 
problems (Sec. |5]l. In Section|6]we outline the weak formulation of Maxwell's equations which is needed 
for the discretization of Maxwell's equations with the finite element method. The basic ideas of the dis- 
cretization are given in Sec. [T] In the discretized version the solution to Maxwell's equations is determined 
in a subspace of the function space which contains the continuous solution. The construction of this sub- 
space and therewith vectorial finite elements is given in Sec. [8] Since the finite size of nano-optical devices 
often has to be taken into account one needs to apply transparent boundary conditions to the computational 
domain. Section|9]explaines our approach. It is based on the implementation of the perfectly matched layer 
(PML) method 121 and allows a certain class of inhomogeneous exterior domains 1161 . The discretization 
scheme of the exterior domain is formulated in the context of the pole condition |T3], which generalizes 
radiation conditions for wave propagation. In the final Section [TO] we apply the finite element method to 
the computation of leaky modes in hollow core photonic crystal fibers. We optimize the fiber design in 
order to minimize radiation losses. 
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Fig. 1 (a) 3D layout of a propagation mode problem. The waveguide structure has an invariance in one spatial dimen- 
sion (we choose the z-direction). (b) 2D cross section of the waveguide structure which is sufficient for computation 
of propagating modes. 

2 Problem Classes Many problems of light propagation can be formulated with time-harmonic Maxwell's 
equations. They can be divided into the following classes: 

• Scattering problems: light scattering and transmission through arbitrary obstacles, e.g. meta materials. 



• Resonance problems: eigenmodes in resonators, e.g. cylindrical cavities, vertical cavity surface emit- 
ting lasers (VCSEL) 

• Propagation mode problems: guided light fields in waveguide structures, e.g. ridge waveguides, 
photonic crystal fibers 

The basic equations for all of these classes are Maxwell's eigenvalue equations, which can be formulated 
as a second order curl curl equation for the electric field: 



The finite element method is used to discretize this differential operator. Corresponding to the above classes 
different mathematical problem formulations arise due to different boundary conditions and unknown 
quantities which we will give in the following 3 sections. Resonance and propagating mode problems 
are eigenvalue problems where eigenmodes of the electric field and corresponding resonance frequencies 
respectively propagating constants have to be determined. In scattering problems an incident field of fixed 
frequency is given and the scattered light field from an arbitrarily shaped object has to be determined. 

The assembling of the finite element system however is very similar for all problem classes since always 
the same operator ([T]i is discretized. After discretizing a scattering problem one has to solve a linear system 
of equations with the assembled finite element matrix and right hand side. After discretizing a propagating 
mode or resonance problem one has to solve an eigenvalue problem for the assembled finite element matrix. 

3 Propagation mode problems The geometry of a waveguide system is invariant in one spatial 
dimension along the fiber, see Fig. Uta). Here we choose the z-direction. Then a propagating mode 
is a solution to the time harmonic Maxwell's equations with frequency lo, which exhibits a harmonic 
dependency in z-direction: 
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Epm(a;, y) and Hpni(a;, y) are the electric and magnetic propagation modes and the parameter is called 
propagation constant. If the permittivity e and permeability yu can be written as: 
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we can split the propagation mode into a transversal and a longitudinal component: 

t^pm(x,y) - ^ E,{x,y) 
Inserting (|2| with ^ and ^ into Maxwell's equations yields: 
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Now we define Ez ~ k^E^ and get the propagation mode problem: 



Problem: 












Find tuples (fcz, E_l, i?z) such that: 
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Eq. dTji is a generalized eigenvalue problem for the propagation constant kz and propagation mode Ep,n (x , y) . 
We get a similar equation for the magnetic field Hpin(a:, y) exchanging e and /i. For the numerical analysis 
of a propagation mode problem (in Sec. [TOl i we furthermore define the effective refractive index rics which 
we will also refer to as eigenvalue: 

k 

ries = 7^ with (10) 

ko 

I. 

ko = T-, 

where Ao is the vacuum wavelength of light. Note that we stated the propagating mode problem ^ on 
IR^- Fig- Clb)- This means that we take the infinite exterior of the waveguide into account, which allows to 
compute bounded as well as leaky modes. Leaky modes thereby model radiation losses from the waveguide 
to the exterior Using the finite element method we therefore have to state transparent boundary conditions 
on F, see Fig. [Hb). We realize these boundary conditions with the PML method which also allows 
inhomogeneous exterior domains. We will explain our implementation of the PML in Sec. |9] 
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Fig. 2 (a) Setup for an electromagnetic resonance problem. The walls of the resonator are perfectly conducting. The 
components of the electric field perpendicular to the normal n of the boundary F therefore have to vanish, (b) Setup 
of a scattering problem. The interior domain Q, contains the scatterer S and is embedded into an infinite exterior R" 
with permittivity e^xt and permeability Hsxt- The incoming electric field Ei„ is entering the interior domain via the 
boundary F and is the source for the electric field inside SI. The scattered light Eout is originated within Q. It is 
therefore strictly outward propagating. 

4 Resonance problems Let us consider a cavity fl with boundary F = d^l as depicted in Fig. 12 a). 
Inside O the electric field has to fulfill Maxwell's equations. On F we have to impose boundary conditions. 
Since in nano-optical applications the coupling of a resonator to the exterior often can not be neglected we 
use transparent boundary conditions. 

Now we are interested in resonance modes and corresponding frequencies in the cavity: 



Problem: 

Find tuples {uj, E) such that: 
Maxwell's equation in interior domain: 




e^^V X ^"^V X E = uj'^E in Q 


(11) 


Maxwell's equation in exterior domain: 




e^^V X ^"^V X 'Eout - oj^'Eout = in M" \ 


(12) 


Boundary condition at F: 




(^^^V X EoMt) X n = (/-t^-^V X E) X n 


(13) 


Bout X n = E X n, 


(14) 


Silver-Miiller radiation condition ("boundary condition at infinity") for exterior field: 




lim r fv X Eo„t(r) x ro - ^itVffffAW^ ^ E„„t(r)) = 

r—*oc y C J 


(15) 


uniformly continuous in each direction ro, 




where r the coordinate vector in M", r its norm and ro = p 





Since all fields E which are gradients of a scalar potential $ lie in the kernel of the curl operator, i.e. 
VxE = 0=>E = V$, the above problem has many solutions to w = 0. For a numerical method it 
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is important to guarantee that the approximated (discrete) fields which correspond to these gradient fields 
also have uj = 0. For the finite element method this means that one has to construct carefully appropri- 
ate ansatz function in order to preserve the mathematical structure of Maxwell's equations in the discrete 
version. We will come to this point in more detail in Sec. [T] 

5 Scattering problems The setup of a scattering problem is depicted in Fig. |2b). The region in space 
occupied with the scatterer S is denoted by fl. For simplicity we assumed here that we have a homoge- 
neous exterior with relative permittivity eext and permeability Hext- In the exterior K" we have an incident 
field Ei„ which enters the interior domain fl across its boundary T = d^l and is scattered. The scattered 
field Eout originates inside and is therefore strictly outgoing. The total field is then E = Ei„ + Eout- 
The scattering problem is formulated as follows: 



Problem: 




Maxwell's equation in interior domain: 




e~^V X ^~^V X E - Lu'^E = in 17 


(16) 


Maxwell's equation in exterior domain: 




e~^W X fi^^W X Eout ~ oj^'Eout = in M" \ f7 


(17) 


Boundary condition at F: 




{fj.-^^ X {Ein + Eout)) X n = {iJ-^V X E) X n 


(18) 


{Ein + Eout) X n = E X n, 


(19) 


where the incoming field E^ has to fulfill Maxwell's equations in a neighborhood of the boundary F. 


Silver-MiiUer radiation condition ("boundary condition at infinity") for exterior field: 




lim r fv X Eout{r) x m - x'^^"""*^""* V x Eout{r)) = 
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uniformly continuous in each direction rg, 




where r the coordinate vector in M", r its norm and ~ j,. 





The Silver Miiller radiation condition guarantees that Eout is a strictly outward radiating solution. For 
inhomogeneous exterior domains a generalization of this radiation condition has to be stated. A deeper 
understanding of the meaning of outward radiating fields is made by the pole condition 1131 which char- 
acterizes these fields by the poles of their Laplace transforms. Using the FEM method fl is taken as 
computational domain. On F transparent boundary conditions have to be stated. 

6 Weak formulation of Maxwell's equations For application of the finite element method we have to 
derive a weak formulation of Maxwell's equations. We multiply ( fTSI ) with a vector valued test function 
^ G V = H{curl, f2) Is) and integrate over the domain O: 

{* • [V X ^-ly X E] - w^e * • E} (fr = , V* S (21) 

where bar denotes complex conjugation. After a partial integration we arrive at the weak formulation of 
Maxwell's equations: 
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Fig. 3 (a) Computational domain and (b) triangulation of pliotonic crystal fiber cross section. 

Find 'E eV = H{curl, il) such that 

j |(V X *) • (^"^V X E) - Lo^e'^- e| rfV = y" * • ¥cPr , V* e T^, (22) 

with 

(/^^^V X E) X n = F given on r (Neumann boundary condition). (23) 

Hence one needs Neumann data on F for the electric field in order to solve the weak problem. Therefore 
one has to construct an operator which maps the Dirichlet data of the electric field onto its Neumann 
values respecting the radiation condition. The construction of such a Dirichlet to Neumann operator will 
be explained in Sec. |9] In order to state the weak formulation we define the following bilinear functionals: 

a(w,v) = / (V X w) • {i-r^V X v) - w^e w • V (fir, (24) 

/(w) = j W-F(fr (25) 

The weak formulation of Maxwell's equations then reads: 



Find V E V = H{curl, fl) such that 




a(w, v) = /(w) , Vw G V. 


(26) 



The above equation is an exact reformulation of Maxwell's equations. 



7 Discretization of Maxwell's equations Now we discretize the weak form of Maxwell's equations. 
The finite element method thereby restricts the space to a finite dimensional subspace Vh C V, dimVh = 
Nh- This means the finite element method does not approximate Maxwell's equations itself but the function 
space in which the solution is determined. The finite dimensional finite element space is an approximation 
to the solution space of the continuous problem. The discretized Maxwell's version (compare to ( |26] )) 
simply reads: 
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Find V e V/i C H{curl, fl) such that 




a(w, v) = /(w) , Vw e V/i. 


(27) 



The subspace Vh and therewith the basis for the approximate solution are constructed as follows. One 
starts with a computational domain il for example the cross section of a photonic crystal fiber ©(a). This 
domain is subdivided into small patches, e.g. triangles or quadrilaterals in 2D and tetrahedrons in 3D, Fig 
Ob). On these patches vectorial ansatz functions cpi are defined whose construction will be explained in 
section [8] These functions are usually polynomials of a fixed order whose support is restricted to one or a 
small number of patches. The set {ipi, ... ,ipNf^} of ansatz functions forms a basis of Vh. The approximate 
solution E/i for the electric field is a superposition of these local ansatz functions: 

N 

1=1 

If we insert this expansion of the electric field into the discrete version of Maxwell's equations dZTl ) for 
V and replace Vw E Vh by Vw e {ipi, . . . , (pN,^} (which is equivalent because this is a basis of Vh) the 
discrete Maxwell's equation reads: 

N 

^ aia{(pj , (fii) = f{ipj) , Vj = 1, . . . , iV (29) 

i=l 

which is a linear system of equations for the unknown coefficients af. 

A-d^f (30) 

ai 

with A^i = a{ipj,ipi), fj ^ f{ipj] 

' aN 

The matrix entries a{Lpj, (pi) arise from computing integrals (l24l i. Since the ansatz functions (pi have a 
small support the stiffness matrix Aji has 0{N) nonzeros out of 0{N'^) entries. The arising matrix is 
therefore sparse. Using special solvers the computational time scales practically linearly with the number 
of unknowns. 

8 Construction of finite elements In this section we want to show how appropriate ansatz functions ipi 
of the finite dimensional subspace Vh are constructed for the finite element method. We restrict ourselves 
to the 2D case and consider triangles as patches which we will denote with K. The global function space 
Vh with dimVh = Nh is separated into local function spaces Vk with dimVK = Nk on each patch K. 
For the union of all patches we have VJK = Q,. 

On each patch we now define Nk basis functions ipi, i = 1, . . . , Nk- Furthermore we have to define 
Nk functionals ipj which are called degrees of freedom. These functionals are constructed such that 

iljj{ipi) ^ Sij (31) 

is satisfied. The meaning of the degrees of freedom can be understood when considering the approximation 
of the solution on a patch, which is a superposition of the local functions (pf. 

Nk /Nk \ 

ai(pi V'i X! "'"^^ ^ - ^^^^ 

1=1 \i=l / 

Hence the degree of freedom i/jj returns the coefficient aj of the basis function (pj. 
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Fig. 4 Set of basis functions Ai of function space P^{K) i33h 

Now we start constructing finite elements. Therefore we have to perform the following steps. Choose 
a patch K, e.g. triangle or quadrilateral. Define a function space Vk on K which has some desired 
properties. E.g. if the field which we are approximating is differentiable then the elements of Vk should 
also be differentiable. Finally define the degrees of freedom and a set of basis functions (pi which span 
Vk. 

We start with linear scalar finite elements of order 1 on triangles. These elements are conform. 
A function in has a first weak derivative. Furthermore the function itself and its weak derivative are 
quadratic Lebesgue integrable lH. Let us consider a triangle with nodes {xi,yi), (x2, 2/2), (2^3, Us)- The 
polynomial function space of order 1 on this triangle is 

{K) = {v = a + bx + cy, a,b,c'ER} . (33) 

This is our local ansatz space with dimP^{K) ~ 3. Now we want to construct the basis functions in 
P^{K) which we will call A; and degrees of freedom tpj. First we define the degrees of freedom. For a 
V S P^{K) we choose: 

ipi{v) := / 6 [{x,y) - {xt,yt)]vdxdy = v{xt,yi), i = 1,2,3 
Jk 

hence the degrees of freedom simply give the value of a function of P^ (K) at the nodes of the patch. Using 
V'j('^i) = ^ij ( I3TI ) we can construct a basis of P^{K). Fig. |4]shows the three basis functions A^ on a unit 
triangle. 

The basis {Ai, A2, A3} is called nodal basis since the degrees of freedom are associated to the nodes 
of the patch. Furthermore each basis function A^ has the value 1 at node i of the patch and the value 
at all other nodes. The global ansatz functions of the finite element space Vh are constructed from these 
local functions. Thereby several patches of the discretized domain can share the same node. On each of 
those patches we have a nodal basis function A; which corresponds to the joint node. All of these basis 
functions then have the same degree of freedom. This means that a global ansatz functions consists of all 
local ansatz functions which have the value 1 at the same node. The global ansatz functions are therefore 
globally continuous. 

Next we want to construct vectorial finite elements of lowest order (but not constant) on the triangle 
which are H{curl) conform. This means that the functions itself and the curl of the functions are quadratic 
Lebesgue integrable 0. As already mentioned it is important that the mathematical structure of the dif- 
ferential operators appearing in Maxwell's equations also holds for their discretized versions, i.e. the 
operators acting on our constructed function spaces of finite dimension. For the continuous operators we 
have the following important property: On simply connected subsets Vl of the following exact sequence 
holds: 



(34) 
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where L^(ri) is the set of functions which are quadratic Lebesgue integrable on H, and H^{n)/M. is the 
set of non constant functions in H^{n). This sequence means: the operator V has an empty kernel on 
(il) /R. The range of V is a subset of H{curl, O) and it is exactly the kernel of V x . The range of V x 
is the whole Hence we construct the local function spaces in such a way that the exact sequence: 

W^/R ^Vh^Sh (35) 
holds, where 

Wh/M. c H\i})/R, (36) 

Vh C H{curl,n), (37) 

Sh C L^in). (38) 

Constructing the linear scalar finite elements on a patch K we already found: 

Wh = {we H\n) : w\k e P\K),yK} . (39) 

This can be seen when evaluating 

ViWh/mK = { ( ) ■ ^ e »} ■ (40) 

We have 

dim {Wh/M.) = dim (V(W,,/M)) = 2 (41) 

hence V has an empty kernel on W/j/R. Furthermore the functions in V(Vr?i/M) lie in the kernel of the 
curl operator Now we want to construct V/j. The exact sequence for the discrete spaces (|35] ) tells us that 
V(VFft,/R) C Vh, which are the constant vectors. Since we wanted to have functions of lowest order but 
not only constant functions in Vh we have to extend it. Again the exact sequence tells us how to make this 
extension: 

VxVh = ShdL^m (42) 

With only constant elements in Vh it follows that 5*^ = {0}. Let us extend Vh in such a way that Sh 
includes at least the constant functions Sh = {s £ L^(ri)/K : s\k & P'^{K)}- The vectorial functions in 
Vh which we include in addition to the constant functions are polynomials of x and y. In the a;-component 
the j/-variable is not allowed to have a degree higher than one and vice versa. Otherwise the curl of the 
vector would not be constant. Therefore we use the following ansatz: 



pi{x) + by 
cx+p2{y) 



£ Vh ^ V X V = c - b, (43) 



where pi, p2 are arbitrary non constant polynomials. Since the exact sequence tells us that only elements 
from V(VF/i/R) lie in the kernel of V x we can deduce pi = P2 = 0: vectors of the form ) \ also 



P2iy) 

lie in the kernel of Vx but are not elements of V(VK/i/M). Since we wanted to extend Vh in a minimal 
way we furthermore choose c = —b. We found: 



Vh = N^{K) = I'' = ( ) + ^ ( ) - ^v^b e m| , dim{N^{K)) = 3. 



(44) 



These elements of lowest order were discovered independently by a number of authors, see [S) . The above 
considerations can be extended to higher order elements as well and were first constructed by Nedelec ID. 
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After finding the local function space we define the degrees of freedom. We associate them to the three 
edges connecting the nodes 1 — ^ 2, 2 ^ 3, 3 ^ 1 of the patch: 



V'l2(v) = 
'/'3l(v) = 



12 



23 



V • rds, 



V • rds. 



V • rds. 



(45) 
(46) 
(47) 



31 



where the integral / is carried out along the edge ij which starts at node i and ends at node j. The quantity 

T is the tangential vector of the edge. The degrees of freedom therefore correspond to the integral of the 
tangential component along the edges of the patches. The basis functions which we construct using OTT l 
are therefore associated to the edges of the patch and the finite elements are called edge elements: 



V'12(</'12) = 1 = 



23 



V'3l(<i2l2) 



= 



31 



ay,i 
ax,2 

ax,3 
fly, 3 



y 

— X 



■ rds. 



■ rds. 



■ rds. 



Carrying out the integrals for all basis functions we can determine all unknown coefficients. The resulting 
basis functions are jS): 



<^12 


= A1VA2 - 


A2VA1 




= A2VA3 - 


A3VA2 


f31 


= A3VA1 - 


A1VA3 



where Ai are the nodal basis functions. The basis {(pi2, 'P23, fsi} is shown in Fig. |5] When we construct 
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Fig. 5 Set of basis functions (pij of function space N^{K), Eq. J44b . 

global basis functions from the local functions Lpij then the local functions of neighboring patches which 
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correspond to the joint edge share the same degree of freedom. Since the degree of freedom is defined 
via the integral of the tangential component this means that the global ansatz functions have a continu- 
ous tangential component. The discretization of the computational domain is performed respecting the 
material boundaries (i.e. all material boundaries lie on edges of the discretization). Therefore the finite 
element solution generically includes the physical property of electric fields having a continuous tangential 
component across boundaries. The normal component can be discontinuous in general. 



9 Transparent boundary conditions In this section we present our realization of transparent boundary 
conditions I1I6J . We use the perfectly matched layer (PML) method ||2| and construct the Dirichlet to 
Neumann operator mentioned in Sec. [3] For simplicity and the sake of a clear presentation we restrict 
ourselves to the 2D case. The ideas carry over to the 3D case of Maxwell's equations. The computational 
domain has to be polygonal and star-shaped. A certain class of inhomogeneities in the exterior is allowed 
which also covers open waveguide structures playing an essential role in integrated optics, see Fig. |6l 
Here a computational domain with a waveguide structure is shown. The Silver-Miiller radiation condition 
(I20I 1 does not hold true for such inhomogeneous domains and has to be generalized. The pole condition 
was introduced in llT3l as a general concept to define radiation conditions for wave propagation problems. 
This also gave new insight into the PML method. The implementation is based on prismatoidal coordinate 



Fig. 6 Prismatoidal coordinate system r;) introduced in the exterior domain. The waveguide structure yields 
solutions analytic in ^-direction. 

systems, shown in Fig. |6] New coordinates ^ and rj are introduced, where ^ denotes a generalized distance 
variable. The central idea is to decompose the exterior domain ilexf' ii^to a finite number of segments 
I2] Each of the segments carries its own coordinate system such that a global distance variable ^ can be 
introduced. Fig. |7] Maxwell's equations are then semi-discretized in a generalized angular variable 77 and 
the PML method is realized via the complex continuation along the ^-direction. The decomposition of the 
exterior is based on straight non-intersecting rays which connect each vertex of the polygonal boundary 
dfl with infinity, shown in Fig. |2l The arising semi-infinite quadrilaterals Q^^'^^ with their xy-coordinate 
system can then be given as the range of reference rectangles Qf'^"* under a local bilinear transformation 




' 1 



FigEl 



(48) 
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Fig. 7 Prismatoidal coordinate system. Each segment Qj is the image of a reference element under a bilinear mapping 
These local mappings are combined to a global mapping B which is continuous in r). 



We then have a transformation of the semi-infinite rectangle Qj^'''"* := [0, oo] x [77^ , 77^+1] onto Q^^'^'' such 
that lines with ^ = const remain parallel under The exact mathematical definitions for the coordinate 
system, radial and normal rays are given in |[T6l . 

Now Maxwell's equations are formulated in the f 77 coordinate system: 

V^^ X pi^-^y^rj X E(^^^) - ^^£,£(5^^) = 0, (49) 

with transformed permittivity and permeability: 

e, = \J\J-^eJ-^ (50) 

where J denotes the Jacobian of the coordinate transformation B^°'' and | J| its determinant. In order to 
formulate the semi-discrete formulation of Maxwell's equation one first has to derive the weak formulation 
of Maxwell's equations. Therefore ( |49l ) is multiplied with a test function w and integrated along F = dVL. 
The rigorous mathematical definitions of the involved function spaces is beyond the scope of this paper 
and can be found in 1161 . The semi discretization with respect to rj of the scattered outgoing electric field 
is performed via the ansatz: 

Nb 

KuM,v) = Y.^outAOi'M (52) 

where the functions . . . , ^'Wi,} form a basis of Sh which is the trace space of the finite element 
space Vh of the interior domain. This means that each ijji can be found in the set of basis functions in Vh 
restricted to the boundary F. Using ansatz (|52] | in the weak formulation of Maxwell's equations gives us 
a linear system of differential equations for the coefficient vector Eq„j(^) ~ {E'^^^l{^), . . . , E'^^^ at^ (C))- 
We will denote it by: 

^(OLeL(O = 0, (53) 

where j4(^)^„( is a differential operator acting on the coefficient vector Eq„((^). The PML method is 
reahzed by replacing C ^ 7? i^il) > 0, 3(7) > 0), which is the complex continuation of the exterior 
solution. Furthermore the unbounded domain fJg^'j'' is replaced with the bounded domain ftp ml ■= 
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{(^,77) G ^ext ■ C ^ [0,p]',il ^ [VminjVmax]}- Becausc of the expected absorbing character of the PML 
we impose zero Dirichlet boundary conditions on the outer boundary ^ = p- The resulting PML system 
then reads 

MlOout^PMLiO=0, (54) 

with Epj^j]^{(,) = £'^„j(7^). Finally the PML system ( l54l i is discretized in the ^ variable with the finite 
element method. 

The complete discretization of the exterior domain may then be interpreted as a FEM discretization on 
quadrilaterals. Their quality depends on the initial choice of the rays. The solution in the exterior is analytic 
in ^ direction. It is therefore advantageous to choose high order finite elements for the ^-discretization. 




(c) 




Fig. 8 First, second and fourth fundamental core modes of HCPCF illustrated in Fig. [S^a) - Parameters: A — 1550 nm, 
r = 300 nm, to = 50 nm, t = 170 nm, 6 cladding rings, A = 589 nm, see Fig. ll U a) for definition of parameters. 

10 Application: Optimization of photonic crystal fiber design In the last section we want to apply 
the finite element method to the computation of propagating modes in hollow core photonic crystal fibers 
(HCPCFs) 15] [I2I m. The results which are presented here are a summary of our work published in 
ifTOl [TTl [Tsl . The computational domain and triangulation of a HCPCF was already shown in Fig. [3 a) 
and (b). The mathematical formulation for this problem type was given in Sec. |3] Here we are interested 
in radiation losses from photonic crystal fibers which means that we compute leaky propagating modes . 
Therefore we have to take the exterior of the fiber into account and apply transparent boundary conditions 
to the computational domain. Fig. [8] shows the first, second, and fourth fundamental leaky core modes. 
We approximated the exterior of the computational domain by a glass cladding of infinite size. This 
approximation is justified if no light which is leaving the microstructured core of a HCPCF is reflected back 
from the outside of the cladding. When applying transparent boundary conditions to the computational 
domain the propagation constant kz © becomes complex fHl . The corresponding leaky mode is therefore 
exponentially damped according to exp(— ^(fc^)^) while propagating along the fiber It can be shown from 
Maxwell's equations that the imaginary part of is proportional to the power flux of the electric field 
across the boundary F of the computational domain ifTsl . This is a quantity one often wants to minimize 
in application. We will quantify the radiation losses by the imaginary part of the effective refractive index 
( [Tol l. First we are interested how accurately we can compute the eigenvalues. Fig. |9]shows the convergence 
of the real and imaginary part of the effective refractive index for an uniform and adaptive grid refinement 
strategy and different finite element degrees ifTDl . For both refinement strategies the real part converges very 
fast the imaginary part is much harder to compute iHl . We see that high order finite elements are needed to 
get an accurate solution for the imaginary part. In an uniform refinement step each triangle is subdivided 
into four smaller ones. This allows to compute the FEM solution more accurately and the relative error 
of the real and imaginary part decreases, see Fig. |9ja), (b). On the other hand the number of unknowns 
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uniform: (a) (b) 





Fig. 9 Relative error of fundamental eigenvalue in dependence on number of unknowns of FEM computation for 
uniform and adaptive refinement strategies and finite element degrees p. Parameters: A = 1550 nm, r = 300 nm, 
w = nm, t = 170 nm, 6 cladding rings, wavelength A = 589 nm. 



and therewith the computational time and memory requirements increase. However if one utilizes the 
sparsity of the linear equation system ( [30b which has to be solved the memory and computational time 
scale linearly with the number of unknowns lU. The eigenvalue obtained from the most accurate FEM 
computation is used as the reference solution for the convergence plots. Since for the finite element method 
convergence is proven mathematically it is reasonable to assume that the FEM solution converges towards 
the exact continuous solution. Fig. |9jc), (d) shows convergence for an adaptive refinement strategy. In 
an adaptive refinement step only a part of the triangles are refined. Since the FEM method works on 
irregular meshes the refinement of only a fraction of all triangles offers no principal difficulties. In order 
to choose the triangles which are refined a so called residuum is evaluated on each triangle ifTSll . This 
residuum quantifies the error of the solution on each triangle and only triangles with the largest residuum 
are refined. According to the definition of the residuum (i.e. the measure for the error) different quantities 
of the electric field converge with a high rate. Therefore the grid can be refined goal-oriented. E.g. if 
one is interested in the field energy one defines a residuum such that this quantity converges with a high 
convergence rate. This refinement strategy was chosen in Fig. |9|c), (d). Compared to uniform refinement 
a certain level of accuracy for the effective refractive index can be achieved with a much smaller number 
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(a) (b) 




r 

Fig. 10 Imaginary part of effective refractive index S5(neff) in dependence on: (a) number of cladding rings, (b) pitch 
A, (c) core surround thickness t, (d) strut thickness w, (e) hole edge radius r. Parameters: A = 1550 nm, r = 300 nm, 
lu = 50 nm, f = 170 nm, 6 cladding rings, wavelength A = 589 nm. 

of unknowns. Another example for the quantity of interest used for goal-oriented grid refinement could be 
the imaginary part of the effective refractive index ifTSlfTTl . 

After we have looked at the convergence we want to optimize the design of a HCPCF in order to 
minimize radiation losses ifTOll . The fiber we are considering has a hollow core corresponding to 19 omitted 
hexagonal cladding cells. It is surrounded by hexagonal cladding rings. The cladding rings form a photonic 
crystal structure which prevents leakage of light from the core to the exterior. With an increasing number 
of cladding rings the radiation losses therefore decrease, see Fig. [TOla) lITOl . We fix the number of cladding 
rings to 6. The free geometrical parameters are the pitch A, hole edge radius r, strut thickness w, and core 



16 J. Pomplun, S. Burger, L. Zschiedrich, and F. Schmidt: Adaptive FEM for simulation of optical nano structures 



(a) (b) 




Fig. 11 (a) geometrical parameters describing HCPCF: pitch A, hole edge radius r, strut thickness w, core surround 
thickness t; (b) detail from a triangulation of HCPCF. Due to the flexibility of triangulations all geometrical features 
of the HCPCF are resolved. 



surround thickness t depicted in Fig. fTTTa) together with the triangulation [TTT b) flOl . Since the finite 
element method works on irregular meshes the modelling of a complicated structure offers no difficulties. 
Fig. [To] shows the radiation losses in dependence on the chosen geometrical parameters keeping all but 
one fixed in each plot flO\. While for the strut thickness w and hole edge radius r we only find one 
minimum for the pitch A and core surround thickness t a large number of local minima can be seen. 
Since with the FEM computation we are able to map the geometrical parameters on the radiation loss we 
can also use an optimization algorithm to find a geometry with minimal attenuation. Therefore we only 
have to choose initial values, e.g. from the one dimensional parameter scans Fig. [TO] We fix the hole 
edge radius to r = 354 nm since it has the weakest effect on the radiation losses. For the starting values 
A = 1550 nm, t ~ 152 nm, w ~ 50 nm optimization yields a minimum value of 5(noff) — 5 • 10~^^^ 
for the imaginary part of the effective refractive index. The corresponding geometrical parameters are 
A = 1597nm,u! = 38nm, i = 151nm||T0l. 
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